fest.f90 Source File
Perform hydrological simulation
Source Code
!! Perform hydrological simulation !|author: <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a> ! license: <a href="http://www.gnu.org/licenses/">GPL</a> ! !``` ! ________ _______ ________ _________ ! |\ _____\|\ ___ \ |\ ____\ |\___ ___\ ! \ \ \__/ \ \ __/| \ \ \___|_ \|___ \ \_| ! \ \ __\ \ \ \_|/__ \ \_____ \ \ \ \ ! \ \ \_| \ \ \_|\ \ \|____|\ \ \ \ \ ! \ \__\ \ \_______\ ____\_\ \ \ \__\ ! \|__| \|_______| |\_________\ \|__| ! \|_________| ! S P A T I A L L Y - D I S T R I B U T E D ! H Y D R O L O G I C A L - M O D E L !``` !### Description ! _FeST_ is a spatially distributed physically based hydrological ! model. _FeST_ is the acronym of “flash–Flood Event–based Spatially ! distributed rainfall–runoff Transformation”. It was initially developed ! by Mancini (1990), as a model oriented to the simulation of ! rainfall-runoff transformation of single flood events. ! Later the _FeST_ model was merged with the soil water balance scheme ! from TOPLATS model (Famiglietti and Wood, 1994), transforming it into a ! continuous model (Montaldo et al., 2007). Then the _FeST_ code was ! redesigned and rewritten from scratch while keeping the basic assumptions ! of the previous release (Rabuffetti et al., 2008). In 2011 the _FeST_ ! was upgraded with a routine to solve the system of water mass and ! energy balance in order to better simulate the actual evapotranspiration ! and interface the model to remotely sensed data (Corbari et al., 2011; ! Corbari & Mancini, 2014). At the same year, 2011, a new module for ! simulating groundwater flux and river-groundwater interaction ! was developed and implemented in the _FeST_ (Ravazzani et al., 2011). ! In 2013 a new version of the code was released built on top of ! the MOSAICO library (Ravazzani, 2013). In 2014 the _FeST_ ! was upgraded with a module for glaciers modelling (Boscarello et al., 2014). ! In 2021 a forest growth component was implemented in the _FeST_ ! (Feki et al., 2021). ! ! !### References ! ! Boscarello, L., Ravazzani, G., Rabuffetti, D., & Mancini, M.. (2014). ! Integrating glaciers raster-based modelling in large catchments ! hydrological balance: the Rhone case study. ! Hydrological processes, 28, 496–508. ! ! ! Corbari, C., & Mancini, M. (2014). Calibration and validation of a distributed ! energy–water balance model using satellite data of land surface temperature ! and ground discharge measurements. Journal of hydrometeorology, 15(1), 376-392. ! ! Corbari, C., Ravazzani, G., & Mancini, M. (2011). A distributed thermodynamic ! model for energy and mass balance computation: FEST–EWB. ! Hydrological Processes, 25(9), 1443-1452. ! ! Famiglietti, J.S., & Wood, E.F. (1994). Multiscale modeling of spatially ! variable water and energy balanceprocess. Water Resour. Res., ! 30, 3061–3078, https://doi.org/10.1029/94WR01498. ! ! Feki, M., Ravazzani, G., Ceppi, A., Pellicone, G., & Caloiero, T.. (2021). ! Integration of forest growth component in the fest-wb distributed ! hydrological model: the bonis catchment case study. Forests, 12(12). ! ! Mancini, M.. (1990). La modellazione distribuita della risposta idrologica: ! effetti della variabilità spaziale e della scala di rappresentazione ! del fenomeno dell'assorbimento. PhD thesis. Politecnico di Milano, ! Istituto di idraulica ! ! Montaldo, N., Ravazzani, G., & Mancini, M.. (2007). On the ! prediction of the toce alpine basin floods with distributed ! hydrologic models. Hydrological processes, 21, 608–621 ! ! Rabuffetti, D., Ravazzani, G., Corbari, C., & Mancini, M.. (2008). ! Verification of operational quantitative discharge forecast (QDF) ! for a regional warning system – the AMPHORE case studies in the ! upper Po river. Natural hazards and earth system sciences, 8, 161–173. ! ! Ravazzani, G.. (2013). Mosaico, a library for raster based hydrological ! applications. Computers & geosciences, 51, 1–6. ! ! Ravazzani, G., Rametta, D., & Mancini, M.. (2011). Macroscopic ! cellular automata for groundwater modelling: a first approach. ! Environmental modelling & software, 26(5), 634–643. ! ! ! !### License ! license: GNU GPL <http://www.gnu.org/licenses/> ! PROGRAM Fest USE IniLib, ONLY: & !Imported derived types: IniList, & !Imported routines: IniOpen, IniClose, & IniReadInt, IniReadReal, IniReadDouble, & IniReadString, SectionIsPresent, & KeyIsPresent USE HydroNetwork, ONLY : & !imported definitions: !NetworkSurfaceFlow, NetworkGroundwater USE Reservoirs, ONLY : & !imported definitions: Reservoir, & !Imported variables: dtReservoir, & !Imported routines: ReservoirSaveStatus USE Diversions, ONLY : & !imported variables: dtDiversion, & !imported routines: DiversionSaveStatus USE SpatialAverage, ONLY : & !imported routines: InitSpatialAverageMeteo, & InitSpatialAverageBalance, & InitSpatialAverageSnow, & InitSpatialAverageGlaciers, & InitSpatialAverageSediment, & InitSpatialAverageCanopy, & InitSpatialAveragePlants, & ComputeSpatialAverageMeteo,& ComputeSpatialAverageBalance, & ComputeSpatialAverageSnow, & ComputeSpatialAverageGlaciers, & ComputeSpatialAverageSediment, & ComputeSpatialAverageCanopy, & ComputeSpatialAveragePlants, & ExportSpatialAverageMeteo, & ExportSpatialAverageBalance, & ExportSpatialAverageSnow, & ExportSpatialAverageGlaciers, & ExportSpatialAverageSediment, & ExportSpatialAverageCanopy, & ExportSpatialAveragePlants, & timeSpatialAverageMeteo, & timeSpatialAverageBalance,& timeSpatialAverageSnow, & timeSpatialAverageIce, & timeSpatialAverageSediment, & timeSpatialAverageCanopy, & timeSpatialAveragePlants USE BasinAverage, ONLY : & !Imported routines: InitBasinAverage, & ReadPointFileBasinAverage, & ExportBasinAverage, & !Imported variables: dtBasinAverage USE RasterExport, ONLY : & !Imported routines: InitRasterExport, & ExportMaps USE Snow, ONLY : & !imported routines: SnowInit, SnowUpdate, & SnowPointInit, & !Imported variables: dtSnow, meltCoefficient, & swe, waterInSnow, snowMelt, & rainfallRate USE Glacier, ONLY : & !imported routines: GlacierInit, GlacierPointInit, & GlacierUpdate, & !imported variables: dtIce, icewe, waterInIce, iceMelt USE SoilBalance, ONLY : & !imported routines InitSoilBalance, SolveSoilBalance, & !imported variables: soilMoisture, soilMoistureRZ, & soilMoistureTZ, soilSatRZ, & soilSatTZ, sEff, rainFlag, & runoff, infilt, percolation, caprise, & soilDepth, wiltingPoint, fieldCapacity, & psdi, ksat, et, pet, dtSoilBalance, & QinSoilSub, soilSat, balanceError, & interstormDuration, deltaSoilMoisture USE Infiltration, ONLY : & !imported variables: infiltrationModel, SCS_CN, & PHILIPEQ, GREEN_AMPT, & sEff, raincum, cuminf USE Groundwater, ONLY : & !Imported routines: GroundwaterInit, & GroundwaterUpdate, & GroundwaterPointInit, & GroundwaterOutputInit, & GroundwaterOutput, & GroundwaterRiverUpdate, & !imported variables: dtGroundwater, & vadoseDepth, & riverGroundwaterInteract, & riverToGroundwater, & groundwaterToRiver USE DischargeRouting, ONLY : & !imported routines: InitDischargeRouting, DischargeRoute, & DischargePointInit, & !imported variables: Qin, Qout, Pin, Pout, Plat, & dtDischargeRouting, & topWidth, waterDepth USE Irrigation, ONLY : & !Imported routines: IrrigationConfig, IrrigationUpdate, & !Imported variables: dtIrrigation, & irrigationFlux USE Sediment, ONLY: & !Imported routines: SedimentInit, InterrillDetachmentRate, & SedimentRouting, & !imported variables: interrillErosion, routeSediment, QoutSed, deltaSed USE DomainProperties, ONLY: & !imported variables: mask, albedoGround, albedo, & albedo_loaded, & !imported routines: DomainInit USE MorphologicalProperties, ONLY: & !imported variables: dem, dem_loaded, & flowDirection, & streamNetwork, & !imported routines: MorphologyInit USE Meteo, ONLY: & !Imported routines: MeteoInit, MeteoRead , & MeteoPointInit, & !Imported variables: dtMeteo USE Precipitation, ONLY: & !imported variables: precipitationRate USE PrecipitationDaily, ONLY: & !imported variables: precipitationRateDaily USE AirTemperature, ONLY: & !imported variables: temperatureAir USE AirTemperatureDailyMean, ONLY: & !imported variables: temperatureAirDailyMean USE AirTemperatureDailyMax, ONLY: & !imported variables: temperatureAirDailyMax USE AirTemperatureDailyMin, ONLY: & !imported variables: temperatureAirDailyMin USE SolarRadiation, ONLY: & !imported variables: radiation, netRadiation USE AirRelativeHumidity, ONLY: & !imported variables: relHumidityAir USE WindFlux, ONLY: & !imported variables: windSpeed USE Plants, ONLY : & !imported routines: PlantsConfig, PlantsGrow, & PlantsParameterUpdate, & !imported variables: dtPlants, lai, fvcover, & gpp, npp, carbonroot, & carbonstem, carbonleaf, dbh, & plantsHeight, density, stemyield, & rsMinLoaded, updatePlantsParameters USE PlantsInterception, ONLY: & !imported variables: dtCanopyInterception, canopyStorage, & canopyPT, & !imported routines: InterceptionInit, Throughfall, & AdjustPT USE DataTypeSizes, ONLY: & !Imported parameters: short, long, float USE LogLib, ONLY: & !Imported routines: LogInit, LogStop, Catch, & !Imported variables: verbose USE GeoLib, ONLY: & !Imported routines: GeoInit USE Chronos, ONLY: & !Imported routines: GetMinute, GetHour, GetMonth, & GetDay, GetDayOfWeek, & !Imported definitions: DateTime, & !Imported variables: timeString, & !Imported operators: ASSIGNMENT (=), & OPERATOR (+), & OPERATOR (-), & OPERATOR (==), & OPERATOR (>), & OPERATOR (<), & OPERATOR (<=) USE GridOperations, ONLY : & !Imported operators: ASSIGNMENT (=) USE GridLib, ONLY : & !imported definition: grid_real, & !imported routines: NewGrid, ExportGrid, & !imported parameters: ESRI_BINARY USE StringManipulation, ONLY : & !imported routines StringCompact, StringTokenize, & StringToLong USE CronSchedule, ONLY : & !Imported types: CronTab, & !Imported routines: CronParseString, & CronIsTime IMPLICIT NONE !configuration TYPE (IniList) :: festIni CHARACTER (len = 300) :: arg !! command line arguments CHARACTER (len = 1000) :: string LOGICAL :: iniFound !!says if configuration file has been passed as command argument !simulation time TYPE (DateTime) :: time !current step time TYPE (DateTime) :: timeStart !start time TYPE (DateTime) :: timeStop !stop time !meteo TYPE (DateTime) :: timeNewMeteo TYPE (DateTime) :: timeOutMeteo INTEGER (KIND = short) :: dtOutMeteo !irrigation TYPE (DateTime) :: timeNewIrrigation TYPE (DateTime) :: timeOutIrrigation INTEGER (KIND = short) :: dtOutIrrigation !snow TYPE (DateTime) :: timeNewSnow TYPE (DateTime) :: timeOutSnow INTEGER (KIND = short) :: dtOutSnow !glacier TYPE(DateTime) :: timeNewIce TYPE(DateTime) :: timeOutIce INTEGER (KIND = short) :: dtOutIce !canopy interception TYPE(DateTime) :: timeNewCanopyInterception TYPE(grid_real) :: raineff ! throughfall[m/s] INTEGER (KIND = short) :: dtOutCanopy TYPE (DateTime) :: timeOutCanopy !forest dynamic TYPE(DateTime) :: timeNewPlants INTEGER (KIND = short) :: dtOutPlants TYPE (DateTime) :: timeOutPlants !soil balance TYPE(DateTime) timeNewSoilBalance !bilancio energetico TYPE (grid_real)Ts !temperatura supeficiale TYPE (grid_real)Rnetta TYPE (grid_real)Xle TYPE (grid_real)Hse TYPE (grid_real)Ge TYPE (grid_real)Ta_prec !temperatura aria istante precedente TYPE (grid_real)Ts_prec !temperatura suolo istante precedente !output TYPE(DateTime) :: timeOutSoilBalance INTEGER :: dtOutSoilBalance ! groundwater INTEGER (KIND = short) :: dtOutGroundwater TYPE(DateTime) :: timeNewGroundwater TYPE(DateTime) :: timeOutGroundwater !sediment INTEGER :: dtCalcSediment TYPE(DateTime) :: timeUpdateSediment INTEGER :: dtOutSediment TYPE(DateTime) :: timeOutSediment !sediment routing INTEGER :: dtCalcSedimentRouting TYPE(DateTime) :: timeUpdateSedimentRouting INTEGER :: dtOutSedimentRouting TYPE(DateTime) :: timeOutSedimentRouting ! TYPE(rete_osservativa) :: xsOutSedimentRouting ! discharge routing TYPE (DateTime) :: timeNewDischargeRouting TYPE(grid_real)Velocita TYPE (grid_real) storageChannelSurf !basin average TYPE (DateTime) :: timeNewBasinAverage INTEGER (KIND = short) :: dtOutBasinAverage !raster export TYPE (DateTime) :: timeNewRasterExport !hotstart CHARACTER (len=1000) :: pathHotstart LOGICAL :: saveLast LOGICAL :: hotstart CHARACTER (LEN = 300) :: cron TYPE (CronTab) :: cronHotstart !general variables INTEGER (KIND = short) :: dtCalc !calculation time step INTEGER (KIND = short) :: dtArray (4) CHARACTER (len = 1000) :: path_out INTEGER :: k,i,j LOGICAL :: fatta_previ = .FALSE. !--------------------------end of declaration---------------------------------- ! splash screen IF ( verbose ) THEN CALL Splash () END IF !------------------------------------------------------------------------------ ! initialize log !------------------------------------------------------------------------------ CALL LogInit !------------------------------------------------------------------------------ ! initialize cartographic engine !------------------------------------------------------------------------------ CALL GeoInit ('GeoLib.ini') !------------------------------------------------------------------------------ ! check command line arguments to read configuration file name !------------------------------------------------------------------------------ i = 1 iniFound = .FALSE. DO WHILE ( .not. (arg == '') ) CALL Getarg ( i, arg ) SELECT CASE (arg) CASE ( '-inifile' ) i = i + 1 CALL Getarg ( i , string ) iniFound = .TRUE. CASE DEFAULT !case unknown END SELECT i = i + 1 END DO !------------------------------------------------------------------------------ ! read name of configuration file if not yet defined !------------------------------------------------------------------------------ IF ( .NOT. iniFound) THEN WRITE(*,*) 'Insert the name of configuration file and press enter' READ(*,*) string END IF CALL Catch ('info', 'Fest', & 'start configuration from file: ', argument = string ) !------------------------------------------------------------------------------ ! read configuration file !------------------------------------------------------------------------------ IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) 'Loading configuration...' END IF string = "./conf-test/fest.ini" CALL IniOpen (string, festIni) !------------------------------------------------------------------------------ ! configure time !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'time configuration' ) IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring time...' END IF !check whether time is present IF ( .NOT. SectionIsPresent ('time', festIni) ) THEN !section is mandatory CALL Catch ('error', 'Fest', 'time section missing in configuration file' ) END IF !start time IF (KeyIsPresent('start', festIni, section = 'time')) THEN timeString = IniReadString ('start', festIni, section = 'time') ELSE CALL Catch ('error', 'Fest', 'start missing in time' ) END IF timeStart = timeString CALL Catch ('info', 'Fest', 'start time: ', argument = timeString ) !stop time IF (KeyIsPresent('stop', festIni, section = 'time')) THEN timeString = IniReadString ('stop', festIni, section = 'time') ELSE CALL Catch ('error', 'Fest', 'stop missing in time' ) END IF timeStop = timeString CALL Catch ('info', 'Fest', 'stop time: ', argument = timeString ) !------------------------------------------------------------------------------ ! result folder !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'result folder configuration' ) IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...setting folder for results...' END IF !check whether section is present IF ( .NOT. SectionIsPresent ('result', festIni) ) THEN !section is mandatory CALL Catch ('error', 'Fest', 'result section missing in configuration file' ) END IF !folder IF (KeyIsPresent('folder', festIni, section = 'result')) THEN path_out = IniReadString ('folder', festIni, section = 'result') ELSE CALL Catch ('error', 'Fest', 'folder missing in result' ) END IF CALL Catch ('info', 'Fest', 'result folder: ', argument = TRIM(path_out) ) !------------------------------------------------------------------------------ ! configure simulation domain !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'simulation domain configuration' ) IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring simulation domain...' END IF !check whether section is present IF ( .NOT. SectionIsPresent ('domain', festIni) ) THEN !section is mandatory CALL Catch ('error', 'Fest', 'domain section missing in configuration file' ) END IF !find configuration file IF (KeyIsPresent('conf-file', festIni, section = 'domain')) THEN string = IniReadString ('conf-file', festIni, section = 'domain') ELSE CALL Catch ('error', 'Fest', 'conf-file missing in domain' ) END IF !read data CALL DomainInit (string) !------------------------------------------------------------------------------ ! morphological properties !------------------------------------------------------------------------------ IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring morphological properties...' END IF !check whether section is present. Optional section IF ( .NOT. SectionIsPresent ('morphology', festIni) ) THEN !section is mandatory CALL Catch ('info', 'Fest', 'morphological properties configuration' ) END IF !find configuration file IF (KeyIsPresent('conf-file', festIni, section = 'morphology')) THEN string = IniReadString ('conf-file', festIni, section = 'morphology') !read data CALL MorphologyInit (string, mask) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in morphology' ) END IF !------------------------------------------------------------------------------ ! configure meteo input !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'meteo input configuration' ) IF ( verbose ) THEN WRITE (*,*) WRITE (*,*) '...configuring meteorological forcings...' END IF !check whether section is present IF ( .NOT. SectionIsPresent ('meteo', festIni) ) THEN !section is mandatory CALL Catch ('error', 'Fest', 'meteo section missing in configuration file') END IF !time step IF (KeyIsPresent('dt', festIni, section = 'meteo')) THEN dtMeteo = IniReadInt ('dt', festIni, section = 'meteo') ELSE CALL Catch ('error', 'Fest', 'dt missing in meteo' ) END IF !find configuration file IF (KeyIsPresent('conf-file', festIni, section = 'meteo')) THEN string = IniReadString ('conf-file', festIni, section = 'meteo') ELSE CALL Catch ('error', 'Fest', 'conf-file missing in meteo' ) END IF !read data CALL MeteoInit (string, timeStart, mask, dem, dem_loaded, albedo_loaded) !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'meteo')) THEN dtOutMeteo = IniReadInt ('dt-out-spatial', festIni, section = 'meteo') IF ( dtOutMeteo > 0 ) THEN timeOutMeteo = timeStart ELSE timeOutMeteo = timeStart + (-1) END IF ELSE CALL Catch ('warning', 'Fest', 'dt-out-spatial missing in meteo' ) dtOutMeteo = - 1 timeOutMeteo = timeStart + (-1) END IF !point data out IF (KeyIsPresent('out-point-file', festIni, section = 'meteo')) THEN string = IniReadString ('out-point-file', festIni, section = 'meteo') CALL MeteoPointInit (string, path_out, timeStart) END IF timeNewMeteo = timeStart !------------------------------------------------------------------------------ ! irrigation !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'irrigation configuration' ) !check whether section is present IF ( SectionIsPresent ('irrigation', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring irrigation...' END IF !time step IF (KeyIsPresent('dt', festIni, section = 'irrigation')) THEN dtIrrigation = IniReadInt ('dt', festIni, section = 'irrigation') ELSE CALL Catch ('error', 'Fest', 'dt missing in irrigation' ) END IF !dt output IF (KeyIsPresent('dt-out', festIni, section = 'irrigation')) THEN dtOutIrrigation = IniReadInt ('dt-out', festIni, section = 'irrigation') IF ( dtOutIrrigation > 0 ) THEN timeOutIrrigation = timeStart ELSE timeOutIrrigation = timeStart + (-1) END IF ELSE CALL Catch ('error', 'Fest', 'dt-out missing in irrigation' ) END IF !configuration file IF ( dtIrrigation > 0 ) THEN IF (KeyIsPresent('conf-file', festIni, section = 'irrigation')) THEN string = IniReadString ('conf-file', festIni, section = 'irrigation') CALL IrrigationConfig (string, path_out, dtOutIrrigation ) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in irrigation' ) END IF timeNewIrrigation = timeStart END IF ELSE dtIrrigation = 0 END IF !------------------------------------------------------------------------------ ! canopy interception !------------------------------------------------------------------------------ !_________________________________________ ! da implementare che se c'è neve non calcolo intercettazione BO! serve? !________________________ CALL Catch ('info', 'Fest', 'canopy interception configuration' ) !raineff is used even when canopy interception is not simulated CALL NewGrid (raineff, mask, 0.) !check whether section is present IF ( SectionIsPresent ('canopy-interception', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring canopy interception...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'canopy-interception')) THEN dtCanopyInterception = IniReadInt ('dt', festIni, section = 'canopy-interception') ELSE CALL Catch ('error', 'Fest', 'dt missing in canopy-interception' ) END IF IF ( dtCanopyInterception > 0 ) THEN !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'canopy-interception')) THEN string = IniReadString ('conf-file', festIni, section = 'canopy-interception') CALL InterceptionInit ( string, mask ) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in canopy-interception' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'canopy-interception')) THEN dtOutCanopy = IniReadInt ('dt-out-spatial', festIni, section = 'canopy-interception') ELSE CALL Catch ('error', 'Fest', 'dt-out-spatial missing in canopy-interception' ) END IF timeNewCanopyInterception = timeStart IF (dtOutCanopy > 0) THEN timeOutCanopy = timeStart ELSE timeOutCanopy = timeStart + (-1) END IF ELSE dtOutCanopy = 0 timeNewCanopyInterception = timeStart + (-1) timeOutCanopy = timeStart + (-1) END IF ELSE dtCanopyInterception = 0 dtOutCanopy = 0 timeNewCanopyInterception = timeStart + (-1) timeOutCanopy = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! Plants dynamic simulation !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'plants configuration' ) !check whether section is present IF ( SectionIsPresent ('plants', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring plants simulation...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'plants')) THEN dtPlants = IniReadInt ('dt', festIni, section = 'plants') ELSE CALL Catch ('error', 'Fest', 'dt missing in plants' ) END IF !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'plants')) THEN string = IniReadString ('conf-file', festIni, section = 'plants') CALL PlantsConfig (string, mask, timeStart, timeStop) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in plants' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'plants')) THEN dtOutPlants = IniReadInt ('dt-out-spatial', festIni, section = 'plants') ELSE CALL Catch ('warning', 'Fest', 'dt-out-spatial missing in plants' ) dtOutPlants = 0 END IF IF (dtOutPlants > 0) THEN timeOutPlants = timeStart END IF IF ( dtPlants > 0 ) THEN timeNewPlants = timeStart ELSE dtOutPlants = 0 timeNewPlants = timeStart + (-1) timeOutPlants = timeStart + (-1) END IF ELSE dtPlants = 0 dtOutPlants = 0 timeNewPlants = timeStart + (-1) timeOutPlants = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! snow accumulation and melting !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'snow configuration' ) !rainfallRate is used even when snow is not simulated CALL NewGrid (rainfallRate, mask, 0.) !check whether section is present IF ( SectionIsPresent ('snow', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring snow accumulation and melting...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'snow')) THEN dtSnow = IniReadInt ('dt', festIni, section = 'snow') ELSE CALL Catch ('error', 'Fest', 'dt missing in snow' ) END IF IF (dtSnow > 0) THEN !check dt IF ( dtSnow /= dtMeteo ) THEN dtSnow=dtMeteo WRITE(*,*)'WARNING: dtSnow \E8 posto uguale a dt_acquisizione_meteo' END IF timeNewSnow = timeStart !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'snow')) THEN string = IniReadString ('conf-file', festIni, section = 'snow') CALL SnowInit (string, mask, timeStart ) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in snow' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'snow')) THEN dtOutSnow = IniReadInt ('dt-out-spatial', festIni, section = 'snow') ELSE CALL Catch ('error', 'Fest', 'dt-out-spatial missing in snow' ) END IF IF ( dtOutSnow > 0 ) THEN timeOutSnow = timeStart END IF !point data out IF (KeyIsPresent('out-point-file', festIni, section = 'snow')) THEN string = IniReadString ('out-point-file', festIni, section = 'snow') CALL SnowPointInit ( string, path_out, timeStart ) END IF ELSE dtOutSnow = 0 timenEWSnow = timeStart + (-1) timeOutSnow = timeStart + (-1) END IF ELSE dtSnow = 0 dtOutSnow = 0 timenEWSnow = timeStart + (-1) timeOutSnow = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! glacier !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'glacer configuration' ) !check whether section is present IF ( SectionIsPresent ('glacier', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring glacier accumulation and ablation...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'glacier')) THEN dtIce = IniReadInt ('dt', festIni, section = 'glacier') ELSE CALL Catch ('error', 'Fest', 'dt missing in glacier' ) END IF IF (dtIce > 0) THEN !check dt IF ( dtIce /= dtMeteo ) THEN dtIce = dtMeteo CALL Catch ('info', 'Fest', 'dtIce set equal to dtMeteo' ) END IF timeNewIce = timeStart !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'glacier')) THEN string = IniReadString ('conf-file', festIni, section = 'glacier') CALL GlacierInit (string, mask, timeStart ) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in glacier' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'glacier')) THEN dtOutIce = IniReadInt ('dt-out-spatial', festIni, section = 'glacier') ELSE CALL Catch ('error', 'Fest', 'dt-out-spatial missing in glacier' ) END IF IF ( dtOutIce > 0 ) THEN timeOutIce = timeStart END IF !point data out IF (KeyIsPresent('out-point-file', festIni, section = 'glacier')) THEN string = IniReadString ('out-point-file', festIni, section = 'glacier') CALL GlacierPointInit ( string, path_out, timeStart ) END IF ELSE dtOutIce = 0 timeNewIce = timeStart + (-1) timeOutIce = timeStart + (-1) END IF ELSE dtIce = 0 dtOutIce = 0 timeNewIce = timeStart + (-1) timeOutIce = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! Soil balance !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'soil balance configuration' ) !check whether section is present IF ( SectionIsPresent ('soil-balance', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring soil water balance...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'soil-balance')) THEN dtSoilBalance = IniReadInt ('dt', festIni, section = 'soil-balance') ELSE CALL Catch ('error', 'Fest', 'dt missing in soil-balance' ) END IF IF ( dtSoilBalance > 0 ) THEN timeNewSoilBalance = timeStart !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'soil-balance')) THEN string = IniReadString ('conf-file', festIni, section = 'soil-balance') CALL InitSoilBalance (string, flowDirection, timeStart) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in soil-balance' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'soil-balance')) THEN dtOutSoilBalance = IniReadInt ('dt-out-spatial', festIni, & section = 'soil-balance') ELSE dtOutSoilBalance = 0 END IF IF (dtOutSoilBalance > 0) THEN timeOutSoilBalance = timeStart END IF ELSE dtOutSoilBalance = 0 timeNewSoilBalance = timeStart + (-1) timeOutSoilBalance = timeStart + (-1) END IF ELSE dtSoilBalance = 0 dtOutSoilBalance = 0 timeNewSoilBalance = timeStart + (-1) timeOutSoilBalance = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! groundwater !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'groundwater configuration' ) !check whether section is present IF ( SectionIsPresent ('groundwater', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring groundwater...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'groundwater')) THEN dtGroundwater = IniReadInt ('dt', festIni, section = 'groundwater') ELSE CALL Catch ('error', 'Fest', 'dt missing in groundwater' ) END IF IF ( dtGroundwater > 0 ) THEN timeNewGroundwater = timeStart !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'groundwater')) THEN string = IniReadString ('conf-file', festIni, section = 'groundwater') CALL GroundwaterInit (string) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in groundwater' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-aquifer', festIni, section = 'groundwater')) THEN dtOutGroundwater = IniReadInt ('dt-out-aquifer', festIni, section = 'groundwater') IF ( dtOutGroundwater > 0 ) THEN timeOutGroundwater = timeStart CALL GroundwaterOutputInit ( path_out) ELSE timeOutGroundwater = timeStart + (-1) END IF ELSE CALL Catch ('error', 'Fest', 'dt-out-aquifer missing in groundwater' ) END IF !point data out IF (KeyIsPresent('out-point-file', festIni, section = 'groundwater')) THEN string = IniReadString ('out-point-file', festIni, section = 'groundwater') CALL GroundwaterPointInit ( string, path_out, timestart ) END IF ELSE dtOutGroundwater = 0 timeNewGroundwater = timeStart + (-1) timeOutGroundwater = timeStart + (-1) END IF ELSE dtGroundwater = 0 dtOutGroundwater = 0 timeNewGroundwater = timeStart + (-1) timeOutGroundwater = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! soil erosion and sediment transport !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'sediment erosion configuration' ) !check whether section is present IF ( SectionIsPresent ('sediment', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring soil erosion and sediment transport...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'sediment')) THEN dtCalcSediment = IniReadInt ('dt', festIni, section = 'sediment') ELSE CALL Catch ('error', 'Fest', 'dt missing in sediment' ) END IF IF ( dtCalcSediment > 0 ) THEN !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'sediment')) THEN string = IniReadString ('conf-file', festIni, section = 'sediment') CALL SedimentInit (string, dtCalcSedimentRouting, string) timeUpdateSediment = timeStart timeUpdateSedimentRouting = timeStart ELSE CALL Catch ('error', 'Fest', 'conf-file missing in sediment' ) END IF !spatial data out dt IF (KeyIsPresent('dt-out-spatial', festIni, section = 'sediment')) THEN dtOutSediment = IniReadInt ('dt-out-spatial', festIni, section = 'sediment') ELSE CALL Catch ('error', 'Fest', 'dt-out-spatial missing in sediment' ) END IF IF (dtOutSediment > 0) THEN timeOutSoilBalance = timeStart END IF timeOutSediment = timeStart ELSE dtOutSediment = 0 timeUpdateSediment = timeStart + (-1) timeOutSediment = timeStart + (-1) END IF ELSE dtCalcSediment = 0 dtOutSediment = 0 timeUpdateSediment = timeStart + (-1) timeOutSediment = timeStart + (-1) END IF !TODO ! IF (routeSediment) THEN ! !set file for output of sediment routing ! fileunit_out_sediment_routing = GetUnit () ! OPEN (fileunit_out_sediment_routing,file=string) ! WRITE(*,*)'nome file anagrafica out punti propagazione sedimenti: ',string(1:LEN_TRIM(string)) ! CALL leggi_anagrafica(xsOutSedimentRouting,fileunit_out_sediment_routing) ! !check dt out sediment routing ! IF(.not.(mod(xsOutSedimentRouting%cadenza_misure,dtCalcSedimentRouting)==0)) THEN ! xsOutSedimentRouting%cadenza_misure=dtCalcSedimentRouting ! WRITE(*,*)'WARNING: xsOutSedimentRouting%cadenza_misure non multiplo di dtCalcSedimentRouting' ! WRITE(*,*)' si impone: xsOutSedimentRouting%cadenza_misure = dtCalcSedimentRouting' ! END IF ! WRITE(*,*)xsOutSedimentRouting%num_staz,' sezioni per risultati propagazione sedimenti' ! DO k=1,xsOutSedimentRouting%num_staz ! xsOutSedimentRouting%staz_oss(k)%misura=0. ! END DO ! CALL righe_colonne(xsOutSedimentRouting,grid_r=QoutSed) ! CLOSE (fileunit_out_sediment_routing) ! ! !write header in output file ! fileunit_out_sediment_routing = GetUnit () ! OPEN (fileunit_out_sediment_routing,file=path_out(1:LEN_TRIM(path_out))//'sediment_routing.fest') ! WRITE(fileunit_out_sediment_routing,'(a)')'FEST-EWB: stima dei sedimentogrammi' ! WRITE(fileunit_out_sediment_routing,'(a)')'anagrafica' ! CALL scrivi_anagrafica(xsOutSedimentRouting,fileunit_out_sediment_routing) ! WRITE(fileunit_out_sediment_routing,'(a)') ! WRITE(fileunit_out_sediment_routing,'(a)') ! WRITE(fileunit_out_sediment_routing,'(a)')'dati' ! WRITE(fileunit_out_sediment_routing,'(a," ",\)')'#' ! DO k=1,xsOutSedimentRouting%num_staz-1 ! WRITE(fileunit_out_sediment_routing,'(a," ",\)')(xsOutSedimentRouting%staz_oss(k)%denominazione) ! END DO ! WRITE(fileunit_out_sediment_routing,'(a)')(xsOutSedimentRouting%staz_oss(xsOutSedimentRouting%num_staz)%denominazione) ! ! !time ! timeOutSedimentRouting = timeStart !END IF ! !------------------------------------------------------------------------------ ! spatial average !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'spatial-average configuration' ) !check whether section is present IF ( SectionIsPresent ('spatial-average', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring spatial average for output...' END IF !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'spatial-average')) THEN string = IniReadString ('conf-file', festIni, section = 'spatial-average') CALL InitSpatialAverageMeteo (string, path_out, temperatureAir, & temperatureAirDailyMean, temperatureAirDailyMax, & temperatureAirDailyMin, precipitationRate, & relHumidityAir, radiation,netRadiation, & windSpeed, precipitationRateDaily, irrigationFlux ) CALL InitSpatialAverageBalance (string, path_out, soilMoisture, & soilMoistureRZ, soilMoistureRZ, & runoff, infilt, percolation, & et, pet, caprise, balanceError) CALL InitSpatialAverageSnow (string, path_out, rainfallRate, swe, & meltCoefficient, waterInSnow, snowMelt ) CALL InitSpatialAverageGlaciers (string, path_out, icewe, waterInIce, iceMelt ) CALL InitSpatialAverageSediment (string, path_out, interrillErosion) CALL InitSpatialAverageCanopy (string, path_out, canopyStorage, raineff, canopyPT ) CALL InitSpatialAveragePlants (string, path_out, lai, gpp, npp, carbonstem, & carbonroot, carbonleaf, fvcover, dbh, plantsHeight, & density, stemyield ) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in spatial-average' ) END IF !initialize meteo output file IF ( dtOutMeteo > 0 ) THEN CALL ComputeSpatialAverageMeteo (dtMeteo, temperatureAir, & temperatureAirDailyMean, & temperatureAirDailyMax, temperatureAirDailyMin, & precipitationRate, relHumidityAir, radiation, & netRadiation, windSpeed, precipitationRateDaily, & irrigationFlux) CALL ExportSpatialAverageMeteo ( init = .TRUE. ) END IF !initialize snow output file IF ( dtOutSnow > 0 ) THEN CALL ComputeSpatialAverageSnow ( dtOutSnow, rainfallRate, & swe, meltCoefficient, waterInSnow, & snowMelt ) CALL ExportSpatialAverageSnow (init = .TRUE.) END IF !initialize glaciers output file IF ( dtOutIce > 0 ) THEN CALL ComputeSpatialAverageGlaciers (dtOutIce, icewe, & waterInIce, iceMelt ) CALL ExportSpatialAverageGlaciers (init = .TRUE.) END IF !initailize vegetation canopy output file IF (dtOutCanopy > 0) THEN CALL ExportSpatialAverageCanopy (init = .TRUE.) END IF !initailize plants output file IF (dtOutPlants > 0) THEN CALL ExportSpatialAveragePlants (init = .TRUE.) END IF !initailize soil balance output file IF (dtOutSoilBalance>0) THEN CALL ExportSpatialAverageBalance (init = .TRUE.) END IF !initailize sediment output file IF (dtCalcSediment > 0) THEN CALL ExportSpatialAverageSediment (init = .TRUE.) END IF END IF !------------------------------------------------------------------------------ ! basin average !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'basin-average configuration' ) !check whether section is present IF ( SectionIsPresent ('basin-average', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring basin average for output...' END IF !output file IF (KeyIsPresent('out-point-file', festIni, section = 'basin-average')) THEN string = IniReadString ('out-point-file', festIni, section = 'basin-average') CALL ReadPointFileBasinAverage ( string ) ELSE CALL Catch ('error', 'Fest', 'out-point-file missing in basin-average' ) END IF !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'basin-average')) THEN string = IniReadString ('conf-file', festIni, section = 'basin-average') CALL InitBasinAverage (string, path_out, temperatureAir, & temperatureAirDailyMean, temperatureAirDailyMax, & temperatureAirDailyMin, precipitationRate, & relHumidityAir, radiation,netRadiation, & windSpeed, precipitationRateDaily, irrigationFlux, & swe, soilMoisture, runoff, infilt, percolation, & et, pet, deltaSoilMoisture, snowMelt, icewe, & iceMelt) !set first export time timeNewBasinAverage = timeStart ELSE CALL Catch ('error', 'Fest', 'conf-file missing in basin-average' ) END IF ELSE timeNewBasinAverage = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! raster export !------------------------------------------------------------------------------ CALL Catch ('info', 'FeST', 'raster-export configuration' ) !check whether section is present IF ( SectionIsPresent ('raster-export', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring raster export for output...' END IF !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'raster-export')) THEN string = IniReadString ('conf-file', festIni, section = 'raster-export') CALL InitRasterExport (string, temperatureAir, precipitationRate, & relHumidityAir, radiation, netRadiation, & windSpeed,swe, soilMoisture, runoff, infilt, & percolation, et, pet ) !set first time timeNewRasterExport = timeStart ELSE CALL Catch ('error', 'Fest', 'conf-file missing in raster-export' ) END IF ELSE timeNewRasterExport = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! discharge routing !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'surface routing configuration' ) !check whether section is present IF ( SectionIsPresent ('discharge-routing', festIni) ) THEN !section is optional IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring discharge routing...' END IF !dt IF (KeyIsPresent('dt', festIni, section = 'discharge-routing')) THEN dtDischargeRouting = IniReadInt ('dt', festIni, section = 'discharge-routing') ELSE CALL Catch ('error', 'Fest', 'dt missing in discharge-routing' ) END IF IF ( dtDischargeRouting > 0 ) THEN timeNewDischargeRouting = timeStart !configuration file IF (KeyIsPresent('conf-file', festIni, section = 'discharge-routing')) THEN string = IniReadString ('conf-file', festIni, section = 'discharge-routing') CALL InitDischargeRouting (fileini = string, time = timeStart, & path_out = path_out, & flowdirection = flowDirection, & domain = mask , & storage = storageChannelSurf, & netRainfall = runoff, & dtrk = dtReservoir ) ELSE CALL Catch ('error', 'Fest', 'conf-file missing in discharge-routing' ) END IF !output data file IF (KeyIsPresent('out-point-file', festIni, section = 'discharge-routing')) THEN string = IniReadString ('out-point-file', festIni, section = 'discharge-routing') CALL DischargePointInit (string, path_out, timeStart) END IF ELSE timeNewDischargeRouting = timeStart + (-1) END IF ELSE dtDischargeRouting = 0 timeNewDischargeRouting = timeStart + (-1) END IF !------------------------------------------------------------------------------ ! hot start !------------------------------------------------------------------------------ CALL Catch ('info', 'Fest', 'save-hot-start configuration' ) !check whether section is present IF ( SectionIsPresent ('save-hot-start', festIni) ) THEN !section is optional hotstart = .TRUE. IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '...configuring hot start...' END IF !time IF (KeyIsPresent('time', festIni, section = 'save-hot-start')) THEN cron = IniReadString ('time', festIni, section = 'save-hot-start') CALL CronParseString (cron, cronHotstart ) !CALL ConfigureHotStart (cron) ELSE CALL Catch ('error', 'Fest', 'time missing in save-hot-start' ) END IF !folder IF (KeyIsPresent('folder', festIni, section = 'save-hot-start')) THEN pathHotstart = IniReadString ('folder', festIni, section = 'save-hot-start') ELSE CALL Catch ('error', 'Fest', 'folder missing in save-hot-start' ) END IF !save last time IF (KeyIsPresent('save-last', festIni, section = 'save-hot-start')) THEN IF ( IniReadInt ('save-last', festIni, section = 'save-hot-start') == 1 ) THEN saveLast = .TRUE. ELSE saveLast = .FALSE. END IF ELSE !set to default saveLast = .FALSE. END IF !IF ( dtHotstart > 0 ) THEN ! ! timeNewHotstart = timeStart ! ! !folder ! IF (KeyIsPresent('folder', festIni, section = 'save-hot-start')) THEN ! pathHotstart = IniReadString ('folder', festIni, section = 'save-hot-start') ! ELSE ! CALL Catch ('error', 'Fest', 'folder missing in save-hot-start' ) ! END IF ! !ELSE ! timeNewHotstart = timeStart + (-1) !END IF !ELSE ! dtHotstart = 0 ! timeNewHotstart = timeStart + (-1) ELSE hotstart = .FALSE. END IF IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) '....configuration completed!' END IF !------------------------------------------------------------------------------ ! set calculation time step !------------------------------------------------------------------------------ dtArray (1) = dtMeteo dtArray (2) = dtSoilBalance dtArray (3) = dtGroundwater dtArray (4) = dtDischargeRouting dtCalc = HUGE ( dtCalc ) DO i = 1, SIZE (dtArray) IF ( dtArray (i) > 0 .AND. dtArray (i) < dtCalc ) THEN dtCalc = dtArray (i) END IF END DO !------------------------------------------------------------------------------ ! start simulation !------------------------------------------------------------------------------ IF ( verbose ) THEN WRITE(*,*) WRITE(*,*) 'starting simulation...' WRITE(*,*) WRITE(*,'(a,I10,a)') 'computation time step: ', dtCalc, ' second' END IF time = timeStart !DO WHILE( .NOT. ( time < timeStop) ) DO WHILE( time <= timeStop ) timeString = time IF ( verbose ) THEN WRITE(*,'(a25)') timeString END IF !----------------------------------------------------------------------------- ! read meteorological forcings !----------------------------------------------------------------------------- IF (time == timeNewMeteo) THEN IF ( verbose ) THEN WRITE(*,*)'read meteorological forcings' END IF CALL MeteoRead (time,dem, albedo) timeNewMeteo = timeNewMeteo + dtMeteo IF ( dtOutMeteo > 0 ) THEN timeSpatialAverageMeteo = time CALL ComputeSpatialAverageMeteo (dtMeteo, temperatureAir, & temperatureAirDailyMean, & temperatureAirDailyMax, temperatureAirDailyMin, & precipitationRate, relHumidityAir, radiation, & netRadiation, windSpeed, precipitationRateDaily, & irrigationFlux) END IF END IF !----------------------------------------------------------------------------- ! update irrigation !----------------------------------------------------------------------------- IF (time == timeNewIrrigation) THEN IF ( verbose ) THEN WRITE(*,*) 'update irrigation' END IF CALL IrrigationUpdate (time, soilSat, Qin, dtOutIrrigation, & timeOutIrrigation) timeNewIrrigation = timeNewIrrigation + dtIrrigation END IF !----------------------------------------------------------------------------- ! snow !----------------------------------------------------------------------------- IF (dtSnow>0) THEN IF (time == timeNewSnow) THEN IF ( verbose ) THEN WRITE(*,*)'snow' END IF CALL SnowUpdate (time, mask) timeNewSnow = timeNewSnow + dtSnow timeSpatialAverageSnow = time CALL ComputeSpatialAverageSnow (dtOutSnow, rainfallRate, swe, & meltCoefficient, waterInSnow, & snowMelt) END IF ELSE rainfallRate = precipitationRate END IF !----------------------------------------------------------------------------- ! glacier !----------------------------------------------------------------------------- IF ( time == timeNewIce ) THEN IF ( verbose ) THEN WRITE(*,*) 'glaciers' END IF CALL GlacierUpdate (time, mask, rainfallRate ) timeNewIce = timeNewIce + dtIce timeSpatialAverageIce = time CALL ComputeSpatialAverageGlaciers (dtOutIce, icewe, & waterInIce, iceMelt) END IF !----------------------------------------------------------------------------- !CANOPY INTERCEPTION !----------------------------------------------------------------------------- IF ( dtCanopyInterception > 0 ) THEN IF ( time == timeNewCanopyInterception) THEN IF ( verbose ) THEN WRITE(*,*) 'canopy interception' END IF CALL Throughfall (rainfallRate , lai, mask, fvcover, raineff) !CALL AdjustPT(fvcover, mask, tp_oss, dtetp) timeNewCanopyInterception = timeNewCanopyInterception + dtCanopyInterception timeSpatialAverageCanopy = time CALL ComputeSpatialAverageCanopy (dtCanopyInterception, canopyStorage, raineff, canopyPT) END IF ELSE !copy rainfallRate into raineff raineff = rainfallRate END IF !------------------------------------------------------------------------------------------------------------- !FOREST DYNAMIC !------------------------------------------------------------------------------------------------------------- !check if plants parameters need update IF ( updatePlantsParameters ) THEN CALL PlantsParameterUpdate ( time ) END IF IF ( dtPlants > 0 ) THEN IF ( time == timeNewPlants) THEN IF ( verbose ) THEN WRITE(*,*) 'plants' END IF ! CALL PlantsGrow (t, radiation, temperatureAir, soilMoisture, fieldCapacity, wiltingPoint, relHumidityAir, co2 % npd (1) % value) timeNewPlants = timeNewPlants + dtPlants timeSpatialAveragePlants = time CALL ComputeSpatialAveragePlants (dtPlants, lai, gpp, npp, carbonstem, carbonroot, carbonleaf, fvcover, dbh, plantsHeight, density, stemyield) END IF END IF !----------------------------------------------------------------------------- ! soil water balance !----------------------------------------------------------------------------- IF ( dtSoilBalance > 0 ) THEN IF ( time == timeNewSoilBalance ) THEN IF ( verbose ) THEN WRITE(*,*)'soil water balance' END IF CALL SolveSoilBalance ( time, raineff, irrigationFlux, & flowDirection, fvcover, vadoseDepth) timeSpatialAverageBalance = time CALL ComputeSpatialAverageBalance (dtSoilBalance, soilMoisture, & soilMoistureRZ, soilMoistureTZ, & runoff, infilt, percolation, & et, pet, caprise, balanceError) timeNewSoilBalance = timeNewSoilBalance + dtSoilBalance END IF ELSE IF ( ALLOCATED (runoff % mat) ) THEN DO i = 1, raineff % idim DO j = 1, raineff % jdim runoff % mat(i,j) = raineff % mat(i,j) !percolation%mat(i,j)=0.0 END DO END DO END IF END IF !-------------------------------------------------------------------------- ! groundwater !-------------------------------------------------------------------------- IF (time == timeNewGroundwater) THEN IF ( verbose ) THEN WRITE(*,*) 'groundwater' END IF !update river groundwater interaction discharge IF (riverGroundwaterInteract ) THEN CALL GroundwaterRiverUpdate ( waterDepth, topWidth ) END IF !update groundwater head CALL GroundwaterUpdate ( time, QinSoilSub, percolation, caprise ) IF ( time == timeOutGroundwater) THEN CALL GroundwaterOutput (time) timeOutGroundwater = timeOutGroundwater + dtOutGroundwater END IF timeNewGroundwater = timeNewGroundwater + dtGroundwater ENDIF !----------------------------------------------------------------------------- ! discharge routing !----------------------------------------------------------------------------- IF ( time == timeNewDischargeRouting ) THEN IF ( verbose ) THEN WRITE(*,*)'discharge routing' END IF CALL DischargeRoute (dt = dtDischargeRouting, time = time, & flowdir = flowDirection, & runoff = runoff, dtrk = dtReservoir, & riverToGroundwater = riverToGroundwater, & groundwaterToriver = groundwaterToRiver, & storage = storageChannelSurf) timeNewDischargeRouting = timeNewDischargeRouting + dtDischargeRouting END IF !----------------------------------------------------------------------------- ! sediment detachment !----------------------------------------------------------------------------- IF (time == timeUpdateSediment) THEN CALL Catch ('info', 'FEST-EWB', 'update sediment erosion and deposition') CALL InterrillDetachmentRate (raineff) timeUpdateSediment = timeUpdateSediment + dtCalcSediment timeSpatialAverageSediment = time CALL ComputeSpatialAverageSediment (dtCalcSediment, interrillErosion) END IF !sediment routing IF (routeSediment) THEN IF (time == timeUpdateSedimentRouting) THEN CALL SedimentRouting (dtCalcSedimentRouting, Qin, Velocita, waterDepth) !update timeUpdateSedimentRouting timeUpdateSedimentRouting = timeUpdateSedimentRouting + dtCalcSedimentRouting END IF END IF !----------------------------------------------------------------------------- ! write spatial average results !----------------------------------------------------------------------------- !output meteo spatial average IF (time == timeOutMeteo) THEN IF ( verbose ) THEN WRITE(*,*) 'export meteo' END IF timeOutMeteo = timeOutMeteo + dtOutMeteo CALL ExportSpatialAverageMeteo () END IF !output canopy interception spatial average IF ( time == timeOutCanopy ) THEN IF ( verbose ) THEN WRITE(*,*) 'export canopy interception' END IF timeOutCanopy = timeOutCanopy + dtOutCanopy CALL ExportSpatialAverageCanopy () END IF !output plants spatial average IF ( time == timeOutPlants ) THEN IF ( verbose ) THEN WRITE(*,*) 'export plants' END IF timeOutPlants = timeOutPlants + dtOutPlants CALL ExportSpatialAveragePlants () END IF !output snow spatial average IF (time==timeOutSnow) THEN IF ( verbose ) THEN WRITE(*,*) 'export snow' END IF timeOutSnow = timeOutSnow + dtOutSnow CALL ExportSpatialAverageSnow () END IF !output glaciers spatial average IF (time == timeOutIce) THEN IF ( verbose ) THEN WRITE(*,*) 'export glaciers' END IF timeOutIce = timeOutIce + dtOutIce CALL ExportSpatialAverageGlaciers () END IF !output soil balance spatial average IF ( time == timeOutSoilBalance ) THEN IF ( verbose ) THEN WRITE(*,*) 'export soil balance' END IF CALL ExportSpatialAverageBalance () timeOutSoilBalance = timeOutSoilBalance + dtOutSoilBalance END IF !!output ragguagli falda !IF ( time == timeOutGroundwater ) THEN ! IF ( verbose ) THEN ! WRITE(*,*) 'export groundwater' ! END IF ! timeOutGroundwater = timeOutGroundwater+dtOutGroundwater ! !CALL esporta_ragguagli_falda(path_out) !END IF !output sediment spatial average IF (time == timeOutSediment) THEN IF ( verbose ) THEN WRITE(*,*) 'export sediment' END IF timeOutSediment = timeOutSediment + dtOutSediment CALL ExportSpatialAverageSediment () END IF !output sediment routing !IF (time == timeOutSedimentRouting) THEN ! xsOutSedimentRouting % istante = time ! CALL celle2stazioni(QoutSed,xsOutSedimentRouting) ! CALL scrivi_dato(xsOutSedimentRouting,fileunit_out_sediment_routing) ! !t ! timeOutSedimentRouting = timeOutSedimentRouting + xsOutSedimentRouting % cadenza_misure ! !END IF !----------------------------------------------------------------------------- ! basin average !----------------------------------------------------------------------------- IF ( time == timeNewBasinAverage ) THEN IF ( verbose ) THEN WRITE (*,*) 'export basin average' END IF CALL ExportBasinAverage (time, temperatureAir, & temperatureAirDailyMean, temperatureAirDailyMax, & temperatureAirDailyMin, precipitationRate, & relHumidityAir, radiation,netRadiation, & windSpeed, precipitationRateDaily, irrigationFlux, & swe, soilMoisture, runoff, infilt, percolation, & et, pet, deltaSoilMoisture, snowMelt, icewe, iceMelt ) timeNewBasinAverage = timeNewBasinAverage + dtBasinAverage END IF !----------------------------------------------------------------------------- ! raster export !----------------------------------------------------------------------------- IF ( time == timeNewRasterExport ) THEN IF ( verbose ) THEN WRITE (*,*) 'update raster maps' END IF CALL ExportMaps (time, dtCalc, temperatureAir, precipitationRate, & relHumidityAir, radiation, netRadiation, & windSpeed, swe, soilMoisture, runoff, & infilt, percolation, et, pet ) timeNewRasterExport = timeNewRasterExport + dtCalc END IF !----------------------------------------------------------------------------- ! save status !----------------------------------------------------------------------------- IF ( hotstart ) THEN !IF ( IsTimeHotStart (time) ) THEN IF ( CronIsTime (time, cronHotstart) ) THEN IF ( verbose ) THEN WRITE (*,*) 'save state variables for hot start' END IF CALL SaveHotStart (time) END IF END IF !----------------------------------------------------------------------------- ! time forward !----------------------------------------------------------------------------- time = time + dtCalc IF ( time > timeStop ) THEN IF ( hotstart ) THEN IF ( saveLast > 0 ) THEN IF ( verbose ) THEN WRITE (*,*) 'save final state variables for hot start' END IF CALL SaveHotStart END IF END IF END IF END DO !final message IF ( verbose ) THEN WRITE(*,*) 'simulation finished' END IF !----------------------------------------- ! terminate log !----------------------------------------- CALL LogStop () CONTAINS !============================================================================== !| Description: ! Save state variables for starting a new simulation from the current condition ! SUBROUTINE SaveHotStart & ! (time) IMPLICIT NONE TYPE(DateTime), OPTIONAL :: time CHARACTER (LEN = 26) :: prefix CHARACTER (LEN = 300) :: file !-------------------------end of declarations---------------------------------- IF (PRESENT(time)) THEN prefix = time prefix = prefix (1:19) // '_' prefix (14:14) = '-' prefix (17:17) = '-' ELSE prefix = ' ' END IF !snow IF ( dtSnow > 0 ) THEN !swe file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'swe' CALL ExportGrid (layer = swe, fileName = file, & fileFormat = ESRI_BINARY) !water in snow file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'water-snow' CALL ExportGrid (layer = waterInSnow, fileName = file, & fileFormat = ESRI_BINARY) END IF !glaciers IF ( dtIce > 0 ) THEN !ice water equivalent file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'icewe' CALL ExportGrid (layer = icewe, fileName = file, & fileFormat = ESRI_BINARY) !water in ice file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'water-ice' CALL ExportGrid (layer = waterInIce, fileName = file, & fileFormat = ESRI_BINARY) END IF !soil balance IF ( dtSoilBalance > 0 ) THEN !soil saturation root zone file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'sat-rz' CALL ExportGrid (layer = soilSatRZ, fileName = file, & fileFormat = ESRI_BINARY) !soil saturation transmission zone file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'sat-tz' CALL ExportGrid (layer = soilSatTZ, fileName = file, & fileFormat = ESRI_BINARY) !precipitation status file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'rainflag' CALL ExportGrid (layer = rainFlag, fileName = file, & fileFormat = ESRI_BINARY) !interstorm duration file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'interstorm' CALL ExportGrid (layer = interstormDuration, fileName = file, & fileFormat = ESRI_BINARY) IF ( infiltrationModel == SCS_CN ) THEN !soil retention file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'seff' CALL ExportGrid (layer = sEff, fileName = file, & fileFormat = ESRI_BINARY) !cumulative precipitation file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'raincum' CALL ExportGrid (layer = raincum, fileName = file, & fileFormat = ESRI_BINARY) END IF IF ( infiltrationModel == PHILIPEQ .OR. & infiltrationModel == GREEN_AMPT ) THEN !cumulative infiltration file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'cuminf' CALL ExportGrid (layer = cuminf, fileName = file, & fileFormat = ESRI_BINARY) END IF END IF ! discharge routing IF ( dtDischargeRouting > 0 ) THEN !Qin file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'Qin' CALL ExportGrid (layer = Qin, fileName = file, & fileFormat = ESRI_BINARY) !Qout file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'Qout' CALL ExportGrid (layer = Qout, fileName = file, & fileFormat = ESRI_BINARY) !Qlat file = pathHotstart (1:LEN_TRIM (pathHotstart)) // & prefix (1:LEN_TRIM(prefix)) // 'Qlat' CALL ExportGrid (layer = Plat, fileName = file, & fileFormat = ESRI_BINARY) ! diversions IF ( dtDiversion > 0 ) THEN IF ( PRESENT (time) ) THEN CALL DiversionSaveStatus ( pathHotstart, time ) ELSE CALL DiversionSaveStatus ( pathHotstart ) END IF END IF ! reservoirs IF ( dtReservoir > 0 ) THEN IF ( PRESENT (time) ) THEN CALL ReservoirSaveStatus ( pathHotstart, time ) ELSE CALL ReservoirSaveStatus ( pathHotstart ) END IF END IF END IF ! !Falda !IF (dtGroundwater > 0) THEN ! DO I=1,AltezzaPiezometrica%NumeroStrati ! ! CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// & ! pref_s(1:LEN_TRIM(pref_s))//'strato'//CHAR(64+I), & ! matrice_real=AltezzaPiezometrica%piezometria(I)) ! ! IF (interazione_falda_fiume ) THEN ! ! CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// & ! pref_s(1:LEN_TRIM(pref_s))//'baseflow_falda', & ! matrice_real=BaseFlow_falda) ! ! CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// & ! pref_s(1:LEN_TRIM(pref_s))//'ricarica_falda', & ! matrice_real=ricarica_falda_fiume) ! ! END IF ! ENDDO !END IF !sediment !IF (dtCalcSediment > 0) THEN ! CALL write_raster_bin(path_hotstart(1:LEN_TRIM(path_hotstart))// & ! pref_s(1:LEN_TRIM(pref_s))//'deltaSed',matrice_real=deltaSed) !END IF ! !IF (dtCalcSedimentRouting > 0) THEN ! !!to do ! !END IF END SUBROUTINE SaveHotStart SUBROUTINE Splash () WRITE (*,*) !WRITE (*,'(a)') ' ______ ______ _____ _______ ' !WRITE (*,'(a)') ' | ____| ____|/ ____|__ __|' !WRITE (*,'(a)') ' | |__ | |__ | (___ | | ' !WRITE (*,'(a)') ' | __| | __| \___ \ | | ' !WRITE (*,'(a)') ' | | | |____ ____) | | | ' !WRITE (*,'(a)') ' |_| |______|_____/ |_| ' !WRITE (*,'(a)') !WRITE (*,'(a)') '-o--o o--o o-o o--O-o ' !WRITE (*,'(a)') ' | | | | ' !WRITE (*,'(a)') ' O-o O-o o-o | ' !WRITE (*,'(a)') ' | | | | ' !WRITE (*,'(a)') ' o o--o o--o o ' !WRITE (*,'(a)') WRITE (*,'(a)') ' ________ _______ ________ _________ ' WRITE (*,'(a)') '|\ _____\|\ ___ \ |\ ____\ |\___ ___\ ' WRITE (*,'(a)') '\ \ \__/ \ \ __/| \ \ \___|_ \|___ \ \_| ' WRITE (*,'(a)') ' \ \ __\ \ \ \_|/__ \ \_____ \ \ \ \ ' WRITE (*,'(a)') ' \ \ \_| \ \ \_|\ \ \|____|\ \ \ \ \ ' WRITE (*,'(a)') ' \ \__\ \ \_______\ ____\_\ \ \ \__\ ' WRITE (*,'(a)') ' \|__| \|_______| |\_________\ \|__| ' WRITE (*,'(a)') ' \|_________| ' WRITE (*,'(a)') ' ' ! (*,'(a)') ' E N E R G Y - W A T E R - B A L A N C E ' WRITE (*,'(a)') ' S P A T I A L L Y - D I S T R I B U T E D ' WRITE (*,'(a)') ' H Y D R O L O G I C A L - M O D E L ' !WRITE (*,'(a)') ' Flood Event Spatially distributed rainfall-runoff Transformation model' WRITE (*,'(a)') RETURN END SUBROUTINE Splash !------------------------------------------------------------------------------------------------------------- !------------------------------------------------------------------------------------------------------------- END PROGRAM Fest